SUBROUTINE ROMBERG(f,a,b,n,integral)
	IMPLICIT NONE
	DOUBLE PRECISION f,a,b,integral,r(8,8),sum
	INTEGER n,i,j
	EXTERNAL ATURAN_TRAPESIUM
	n = n/2
	if (n.gt.8) then
		n = 8
	end if
	
	r(1,1)=0.5*(b-a)*(f(a)+f(b))
	do i=2,n
		call ATURAN_TRAPESIUM(f,a,b,2**i-1,r(i,1));
	end do
	
	do i = 2,n
		do j = 2,i
			r(i,j) = r(i,j-1)+(r(i,j-1)-r(i-1,j-1))/(4**(j-1)-1)
		end do
	end do
	integral = r(n,n)
END

SUBROUTINE ATURAN_TRAPESIUM(fungsi, a, b, n, I)
	
	IMPLICIT NONE
	
	DOUBLE PRECISION h, x, sigma, fungsi, a, b, I
	INTEGER r, n
	
	h = (b - a) / n
	x = a
	
	I = fungsi(a) + fungsi(b)
	sigma = 0
	r = 1
	
	do r=1,(n-1)
		x = x + h
		sigma = sigma + 2*fungsi(x)
	end do
	
	I = (I + sigma)*h/2
END

SUBROUTINE RUNGE_KUTTA_ORDE_3(f,x0, y0, b, h,output)
	IMPLICIT NONE
	DOUBLE PRECISION f,x0, y0, b, h, output, x, y, k1, k2, k3
	INTEGER r, n
	
	n = (b - x0)/h
	y = y0
	x = x0
	do r=1,n
		k1 = h*f(x,y)
		k2 = h*f(x+h/2,y+k1/2)
		k3 = h*f(x+h, y-k1+2*k2)
		y = y + (k1+4*k2+k3)/6
		x = x+h
	end do
	output = y;
END

SUBROUTINE RUNGE_KUTTA_ORDE_4(f,x0, y0, b, h,output)
	IMPLICIT NONE
	DOUBLE PRECISION f,x0, y0, b, h, output, x, y, k1, k2, k3, k4
	INTEGER r, n
	
	n = (b - x0)/h
	y = y0
	x = x0
	do r=1,n
		k1 = h*f(x,y)
		k2 = h*f(x+h/2,y+k1/2)
		k3 = h*f(x+h/2, y+k2/2)
		k4 = h*f(x+h, y+k3)
		y = y + (k1+2*k2+2*k3+ k4)/6
		x = x+h
	end do
	output = y;
END

SUBROUTINE ADAM_BASHFORTH_MOULTON(f,x0,y0,b,h,output)
	IMPLICIT NONE
	DOUBLE PRECISION x0,y0,b,h,output,x,y,y1,y2,y3,f
	INTEGER r,n
	EXTERNAL RUNGE_KUTTA_ORDE_3
	
	n = (b-x0)/h
	call RUNGE_KUTTA_ORDE_3(f,x0, y0, x0+h, h,y1)
	call RUNGE_KUTTA_ORDE_3(f,x0, y0, x0+2*h, h,y2)
	call RUNGE_KUTTA_ORDE_3(f,x0, y0, x0+3*h, h,y3)
	x = x0 + 3*h
	do r=4,n
		y = y3 + h/24*(-9*f(x-3,y0) + 37*f(x-2*h,y1) - 59*f(x-h,y2) + 55*f(x,y3))
		y = y3 + h/24*(f(x-2*h,y1) - 5*f(x-h,y2) + 19*f(x,y3) + 9*f(x+h,y))
		y0 = y1
		y1 = y2
		y2 = y3
		y3 = y
		x = x+h
	end do
	output = y
END